Update 28 Jan 2017

Summary of findings:

  1. % load and total mass are too collinear to include in the same model.
  2. Treatment and % load are too collinear to include in the same model.
  3. For the remainder of the modeling, I’ll use only % load

Differences between model with Treatment vs. model with % load

wingbeat freq = size Met rate = mass carrying amplitude = percent loading

Susie’s (hypotheses): lm( frequency ~ ITspan + trt order)

lm( metabolism ~ total mass)

lm( stroke amplitude ~ %load + trt) (those are probably the same thing)

** Hypotheses: 1. Wingbeat freq moves around for a number of reasons (tired, cold, added mass) 2. Stroke amplitude rescues wingbeat frequency 3. As you get really close to your max, adding more weight doesn’t seem to cost as much.


Analyses of respirometry data:

  1. How does carrying a load during flight affect bumblebees’?
    1. How does load affect respiratory rate?
    2. How does load affect wingbeat frequency?
    3. How does load affect stroke amplitude?
    4. How does load affect wing velocity?
  2. Which flight measurements are most closely associated with respiratory rate in bees – wingbeat frequency, stroke amplitude, or wing velocity?

Setup

Install required packages and read in data Define custom function for evaluating VIF with multilevel models

ipak <- function(pkg){
     new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
     if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
     sapply(pkg, require, character.only = TRUE)
}

packages <- c('ggplot2', 'car', 'lme4', 'gsheet', "MASS")
ipak(packages)
## ggplot2     car    lme4  gsheet    MASS 
##    TRUE    TRUE    TRUE    TRUE    TRUE
theme_set(theme_bw())

# read in data -- google sheet called "Bumble mumble grumble"
bdta <- gsheet2tbl("https://docs.google.com/spreadsheets/d/1GUUoFq41Ep3sJNFXiyMcp7fZhYmQCzTcLnhVXdr4WRo/edit?usp=sharing")

Function for calculating vif for lmer

(used later)

vif.mer <- function (fit) {
    ## adapted from rms::vif
    v <- vcov(fit)
    nam <- names(fixef(fit))

    ## exclude intercepts
    ns <- sum(1 * (nam == "Intercept" | nam == "(Intercept)"))
    if (ns > 0) {
        v <- v[-(1:ns), -(1:ns), drop = FALSE]
        nam <- nam[-(1:ns)]
    }
    
    d <- diag(v)^0.5
    v <- diag(solve(v/(d %o% d)))
    names(v) <- nam
    v
}

Data Overview

# calculate percent load
bdta$percLoad <- with(bdta, (Mstarved + load) / Mstarved)

# calculate total mass
bdta$totalMass <- with(bdta, Mstarved + load)

summary(bdta)
##     Bee.ID          Treatment.order  Treatment            Mstarved     
##  Length:60          Min.   :1.0     Length:60          Min.   :0.0767  
##  Class :character   1st Qu.:1.0     Class :character   1st Qu.:0.1213  
##  Mode  :character   Median :1.5     Mode  :character   Median :0.1378  
##                     Mean   :1.5                        Mean   :0.1411  
##                     3rd Qu.:2.0                        3rd Qu.:0.1670  
##                     Max.   :2.0                        Max.   :0.1909  
##       load            IT.Span      single.wing.area av.resp..CO2.mL.hr.
##  Min.   :0.01650   Min.   :3.860   Min.   :20.03    Min.   : 5.507     
##  1st Qu.:0.04713   1st Qu.:4.400   1st Qu.:27.93    1st Qu.: 8.935     
##  Median :0.07840   Median :4.620   Median :30.67    Median :10.666     
##  Mean   :0.08332   Mean   :4.594   Mean   :30.70    Mean   :10.600     
##  3rd Qu.:0.11800   3rd Qu.:4.870   3rd Qu.:33.75    3rd Qu.:12.297     
##  Max.   :0.17460   Max.   :5.180   Max.   :39.05    Max.   :16.475     
##     wbf.aud.     stroke.amplitude wing.velocity      percLoad    
##  Min.   :163.9   Min.   : 96.88   Min.   :1.208   Min.   :1.143  
##  1st Qu.:173.0   1st Qu.:114.16   1st Qu.:1.351   1st Qu.:1.321  
##  Median :179.1   Median :125.96   Median :1.427   Median :1.576  
##  Mean   :178.9   Mean   :123.31   Mean   :1.465   Mean   :1.600  
##  3rd Qu.:185.3   3rd Qu.:133.61   3rd Qu.:1.588   3rd Qu.:1.863  
##  Max.   :194.7   Max.   :144.78   Max.   :1.843   Max.   :2.180  
##    totalMass     
##  Min.   :0.1062  
##  1st Qu.:0.1808  
##  Median :0.2191  
##  Mean   :0.2245  
##  3rd Qu.:0.2646  
##  Max.   :0.3550
ggplot(bdta, aes(x = percLoad, fill = Treatment)) + 
     geom_histogram()  
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# seems like percent load and trt are very related
# visualize % load vs total mass
ggplot(bdta, aes(x = percLoad, y = Mstarved + load)) + 
     geom_point()

# seem to be strongly associated
# see number of observations per bee
table(bdta$Bee.ID) # each bee has two observations
## 
## E01 E03 E05 E09 E10 E11 E12 E16 E20 E24 E26 E28 E30 E31 E32 E33 E35 E36 
##   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2 
## E37 E39 E41 E42 E44 E45 E49 E50 E53 E54 E56 E59 
##   2   2   2   2   2   2   2   2   2   2   2   2
# get number of unque bees
length(table(bdta$Bee.ID))
## [1] 30
# get treatment orders
trtOrders_loadedSecond <- sapply(X = unique(bdta$Bee.ID), FUN = function(x){
     tmp = bdta[bdta$Bee.ID == x, c("Treatment.order", "Treatment")]
     if("2_L" %in% paste(tmp[,1], tmp[,2], sep = "_")){
          loadedSecond = TRUE
     }
     else loadedSecond = FALSE
     return(loadedSecond)
})

table(trtOrders_loadedSecond)
## trtOrders_loadedSecond
## FALSE 
##    30
# visualize scatterplot
car::scatterplotMatrix(bdta[, 4:13])

Variance Inflation Factor and PCA

# VIF  is high among the bee size predictors
# note: this function is almost the same as car::vif
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment + Treatment.order + Mstarved + IT.Span + totalMass +  single.wing.area + percLoad +  (1|Bee.ID), data = bdta))
##      TreatmentUL  Treatment.order         Mstarved          IT.Span 
##        16.591645         1.270430        23.112962         6.658432 
##        totalMass single.wing.area         percLoad 
##        32.931859        10.903672        38.415618
car::scatterplotMatrix(bdta[, c("Mstarved", "IT.Span", "single.wing.area")])

# principle components
aa = prcomp(bdta[, c("Mstarved", "IT.Span", "single.wing.area")], center = TRUE, scale = TRUE)
summary(aa) # 1st pc explains ~95% of the variance in the three measurements of size
## Importance of components:
##                           PC1     PC2     PC3
## Standard deviation     1.6852 0.33843 0.21323
## Proportion of Variance 0.9467 0.03818 0.01516
## Cumulative Proportion  0.9467 0.98484 1.00000
biplot(aa) # shows that all three size measurement are correlated

# note, I changed the signs of the predictions so that higher PC1 values 
# correspond to bigger bees
p1 = -predict(aa)[,1] 

# add PC1 scores to dataset
bdta$size_pca1 = p1

# show scatterplot matrix to see correlations among size predictors
car::scatterplotMatrix(bdta[, c("Mstarved", "IT.Span", "single.wing.area", "size_pca1", "percLoad")])

# check VIF one more time
# VIF  is high among the bee size predictors
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment + Treatment.order + size_pca1 + percLoad + totalMass +  (1|Bee.ID), data = bdta))
##     TreatmentUL Treatment.order       size_pca1        percLoad 
##       16.643643        1.258990        7.448179       30.631628 
##       totalMass 
##       26.790114
# remove treatment 
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment.order + size_pca1 + percLoad + totalMass + (1|Bee.ID), data = bdta))
## Treatment.order       size_pca1        percLoad       totalMass 
##        1.014834        7.698050       20.910198       25.813622
# looks like we cannot have total mass and percentLoad in the same model, according to VIF

# remove total mass 
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment.order + size_pca1 + percLoad +  (1|Bee.ID), data = bdta))
## Treatment.order       size_pca1        percLoad 
##        1.003811        1.004814        1.008624

Multilevel models

# a few house-keeping issues to make models easier to read

# convert trt order to factor
bdta$Treatment.order <- as.factor(as.character(bdta$Treatment.order))

# change reference level to unloaded
bdta$Treatment <- factor(bdta$Treatment, levels = c("UL", "L"))

Resp Rate

# make a full model with all two-way interactions
m1 <- lmer(av.resp..CO2.mL.hr. ~  (size_pca1 +  Treatment.order + percLoad)^2 + (1|Bee.ID), data = bdta)

summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## av.resp..CO2.mL.hr. ~ (size_pca1 + Treatment.order + percLoad)^2 +  
##     (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 196.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5230 -0.4541 -0.0409  0.4813  2.2994 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.2498   1.1179  
##  Residual             0.7711   0.8781  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                            Estimate Std. Error t value
## (Intercept)                  3.0096     1.3291   2.264
## size_pca1                    0.6926     0.4455   1.555
## Treatment.order2             1.9553     2.0573   0.950
## percLoad                     4.8854     0.8030   6.084
## size_pca1:Treatment.order2  -0.0830     0.1451  -0.572
## size_pca1:percLoad           0.3533     0.2491   1.418
## Treatment.order2:percLoad   -1.4903     1.2730  -1.171
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 percLd s_1:T. sz_1:L
## size_pca1   -0.223                                   
## Trtmnt.rdr2 -0.860  0.246                            
## percLoad    -0.980  0.204  0.866                     
## sz_pc1:Tr.2  0.297 -0.379 -0.273 -0.297              
## sz_pc1:prcL  0.148 -0.933 -0.184 -0.125  0.228       
## Trtmnt.r2:L  0.853 -0.233 -0.994 -0.868  0.271  0.169
### LRT's for interactions
m2.0 <- update(m1, .~. - Treatment.order:percLoad)

anova(m1, m2.0) # likelihood ratio test for interaction of treatment order
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.0: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad + 
## m2.0:     (1 | Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## m1: av.resp..CO2.mL.hr. ~ (size_pca1 + Treatment.order + percLoad)^2 + 
## m1:     (1 | Bee.ID)
##      Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.0  8 208.02 224.78 -96.011   192.02                        
## m1    9 208.55 227.40 -95.275   190.55 1.472      1      0.225
summary(m2.0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +  
##     (1 | Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
##    Data: bdta
## 
## REML criterion at convergence: 199.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.61501 -0.46662 -0.04214  0.53248  2.31063 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.3003   1.1403  
##  Residual             0.7595   0.8715  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                            Estimate Std. Error t value
## (Intercept)                 4.33950    0.69115   6.279
## size_pca1                   0.57136    0.43113   1.325
## Treatment.order2           -0.43829    0.22783  -1.924
## percLoad                    4.06743    0.39578  10.277
## size_pca1:Treatment.order2 -0.03698    0.13862  -0.267
## size_pca1:percLoad          0.40235    0.24379   1.650
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 percLd s_1:T.
## size_pca1   -0.048                            
## Trtmnt.rdr2 -0.219  0.137                     
## percLoad    -0.924  0.003  0.055              
## sz_pc1:Tr.2  0.131 -0.337 -0.036 -0.131       
## sz_pc1:prcL  0.007 -0.930 -0.144  0.045  0.192
m2.1 <- update(m2.0,.~. - size_pca1:Treatment.order)
anova(m2.0, m2.1)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.1: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad + 
## m2.1:     (1 | Bee.ID) + size_pca1:percLoad
## m2.0: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad + 
## m2.0:     (1 | Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2.1  7 206.10 220.76 -96.051   192.10                         
## m2.0  8 208.02 224.78 -96.011   192.02 0.0814      1     0.7754
summary(m2.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +  
##     (1 | Bee.ID) + size_pca1:percLoad
##    Data: bdta
## 
## REML criterion at convergence: 197.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.68107 -0.47084 -0.05908  0.51796  2.31287 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.3135   1.1461  
##  Residual             0.7334   0.8564  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)          4.3659     0.6752   6.466
## size_pca1            0.5329     0.4000   1.332
## Treatment.order2    -0.4405     0.2237  -1.969
## percLoad             4.0521     0.3858  10.502
## size_pca1:percLoad   0.4146     0.2352   1.763
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 percLd
## size_pca1   -0.004                     
## Trtmnt.rdr2 -0.216  0.132              
## percLoad    -0.921 -0.044  0.051       
## sz_pc1:prcL -0.019 -0.935 -0.140  0.072
m2.2 <- update(m2.1, .~. - size_pca1:percLoad)
anova(m2.1, m2.2) ## drop all interactions
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.2: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad + 
## m2.2:     (1 | Bee.ID)
## m2.1: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad + 
## m2.1:     (1 | Bee.ID) + size_pca1:percLoad
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m2.2  6 207.38 219.94 -97.688   195.38                           
## m2.1  7 206.10 220.76 -96.051   192.10 3.2723      1    0.07046 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m2.2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +  
##     (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 199.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.44974 -0.62601 -0.00733  0.59569  2.08255 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.296    1.138   
##  Residual             0.785    0.886   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)        4.3843     0.6950   6.309
## size_pca1          1.1922     0.1423   8.378
## Treatment.order2  -0.3854     0.2292  -1.681
## percLoad           4.0055     0.3976  10.074
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2
## size_pca1   -0.064              
## Trtmnt.rdr2 -0.221  0.004       
## percLoad    -0.925  0.069  0.062
# renaming model to simplify later typing
m2 <- m2.2

##### LRTs for main effects
## Treatment Order
m3 <- update(m2, .~. - Treatment.order)
anova(m2, m3, test = "Chi") # drop trt order (different than model with treatment)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m3: av.resp..CO2.mL.hr. ~ size_pca1 + percLoad + (1 | Bee.ID)
## m2: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad + 
## m2:     (1 | Bee.ID)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m3  5 208.25 218.73 -99.127   198.25                           
## m2  6 207.38 219.94 -97.688   195.38 2.8789      1    0.08975 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT for size
m4 <- update(m3, .~. - size_pca1)
anova(m3, m4) # keep size
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m4: av.resp..CO2.mL.hr. ~ percLoad + (1 | Bee.ID)
## m3: av.resp..CO2.mL.hr. ~ size_pca1 + percLoad + (1 | Bee.ID)
##    Df    AIC    BIC   logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m4  4 244.02 252.40 -118.011   236.02                             
## m3  5 208.25 218.73  -99.127   198.25 37.767      1   7.97e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT for Treatment (load)
m5 <- update(m3, .~. - percLoad)
anova(m3, m5) # keep percent load
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m5: av.resp..CO2.mL.hr. ~ size_pca1 + (1 | Bee.ID)
## m3: av.resp..CO2.mL.hr. ~ size_pca1 + percLoad + (1 | Bee.ID)
##    Df    AIC    BIC   logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m5  4 255.27 263.64 -123.634   247.27                             
## m3  5 208.25 218.73  -99.127   198.25 49.013      1  2.543e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# summarize final model for paper
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pca1 + percLoad + (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 201.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.33113 -0.59393  0.02882  0.58497  2.24345 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.2695   1.1267  
##  Residual             0.8349   0.9137  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   4.1216     0.6955   5.926
## size_pca1     1.1932     0.1423   8.388
## percLoad      4.0493     0.4087   9.907
## 
## Correlation of Fixed Effects:
##           (Intr) sz_pc1
## size_pca1 -0.067       
## percLoad  -0.940  0.071
# write output
summary(m3)$coefficients  
##             Estimate Std. Error  t value
## (Intercept) 4.121569  0.6955406 5.925706
## size_pca1   1.193243  0.1422611 8.387692
## percLoad    4.049269  0.4087169 9.907270
write.csv( summary(m3)$coefficients, file = "RespCoefs_percLoad.csv" )


# diagnostics
# qq plot
qqnorm(resid(m3), main = "")
qqline(resid(m3)) # good

# residual plot
plot(fitted(m3), resid(m3), xlab = "fitted", ylab = "residuals")
abline(0,0)

# QQPlot for group-level effects
qqnorm(ranef(m3)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m3)$Bee.ID[[1]]) # looks good

Wingbeat Freq

# make a full model with all two-way interactions
m1 <- lmer(wbf.aud. ~  (size_pca1 +  Treatment.order + percLoad)^2 + (1|Bee.ID), data = bdta)

summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wbf.aud. ~ (size_pca1 + Treatment.order + percLoad)^2 + (1 |  
##     Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 349.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.80184 -0.60216 -0.06242  0.61714  1.64271 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 10.18    3.190   
##  Residual             20.62    4.541   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                            Estimate Std. Error t value
## (Intercept)                173.7585     5.7634  30.149
## size_pca1                   -4.8135     2.1942  -2.194
## Treatment.order2            -5.8644     8.3595  -0.702
## percLoad                     4.3885     3.4953   1.256
## size_pca1:Treatment.order2   0.9035     0.7397   1.221
## size_pca1:percLoad           0.9858     1.2643   0.780
## Treatment.order2:percLoad    1.4493     5.1496   0.281
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 percLd s_1:T. sz_1:L
## size_pca1   -0.198                                   
## Trtmnt.rdr2 -0.818  0.226                            
## percLoad    -0.984  0.172  0.814                     
## sz_pc1:Tr.2  0.254 -0.376 -0.223 -0.252              
## sz_pc1:prcL  0.127 -0.958 -0.171 -0.098  0.216       
## Trtmnt.r2:L  0.807 -0.209 -0.990 -0.818  0.221  0.152
### LRT's for interactions
m2.0 <- update(m1, .~. - size_pca1:percLoad )

anova(m1, m2.0)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.0: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID) + 
## m2.0:     size_pca1:Treatment.order + Treatment.order:percLoad
## m1: wbf.aud. ~ (size_pca1 + Treatment.order + percLoad)^2 + (1 | 
## m1:     Bee.ID)
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2.0  8 381.69 398.44 -182.84   365.69                         
## m1    9 382.99 401.83 -182.49   364.99 0.7005      1     0.4026
summary(m2.0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID) +  
##     size_pca1:Treatment.order + Treatment.order:percLoad
##    Data: bdta
## 
## REML criterion at convergence: 352.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.72938 -0.54168  0.01881  0.54877  1.67034 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 10.08    3.175   
##  Residual             20.49    4.526   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                            Estimate Std. Error t value
## (Intercept)                173.1944     5.6964  30.404
## size_pca1                   -3.1748     0.6277  -5.058
## Treatment.order2            -4.7581     8.2056  -0.580
## percLoad                     4.6533     3.4658   1.343
## size_pca1:Treatment.order2   0.7793     0.7199   1.082
## Treatment.order2:percLoad    0.8413     5.0701   0.166
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 percLd s_1:T.
## size_pca1   -0.268                            
## Trtmnt.rdr2 -0.815  0.221                     
## percLoad    -0.984  0.272  0.813              
## sz_pc1:Tr.2  0.234 -0.606 -0.194 -0.238       
## Trtmnt.r2:L  0.803 -0.222 -0.990 -0.816  0.194
m2.1 <- update(m2.0,.~. - size_pca1:Treatment.order)
anova(m2.0, m2.1)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.1: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID) + 
## m2.1:     Treatment.order:percLoad
## m2.0: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID) + 
## m2.0:     size_pca1:Treatment.order + Treatment.order:percLoad
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2.1  7 380.97 395.64 -183.49   366.97                         
## m2.0  8 381.69 398.44 -182.84   365.69 1.2893      1     0.2562
summary(m2.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID) +  
##     Treatment.order:percLoad
##    Data: bdta
## 
## REML criterion at convergence: 354.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.60845 -0.62485  0.06463  0.64370  1.69662 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 10.35    3.218   
##  Residual             20.39    4.515   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                           Estimate Std. Error t value
## (Intercept)               171.6958     5.5454  30.962
## size_pca1                  -2.7623     0.5021  -5.501
## Treatment.order2           -2.9560     8.0745  -0.366
## percLoad                    5.5797     3.3706   1.655
## Treatment.order2:percLoad  -0.2770     4.9889  -0.056
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 percLd
## size_pca1   -0.163                     
## Trtmnt.rdr2 -0.808  0.133              
## percLoad    -0.983  0.165  0.806       
## Trtmnt.r2:L  0.796 -0.133 -0.989 -0.809
m2.2 <- update(m2.1, .~. - Treatment.order:percLoad)
anova(m2.1, m2.2) ## drop all interactions
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.2: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID)
## m2.1: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID) + 
## m2.1:     Treatment.order:percLoad
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2.2  6 378.98 391.54 -183.49   366.98                         
## m2.1  7 380.97 395.64 -183.49   366.97 0.0018      1     0.9666
summary(m2.2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 360
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.64247 -0.64283  0.06575  0.64332  1.72090 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept)  9.874   3.142   
##  Residual             20.221   4.497   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)      171.9557     3.3390   51.50
## size_pca1         -2.7663     0.4908   -5.64
## Treatment.order2  -3.3999     1.1632   -2.92
## percLoad           5.4191     1.9692    2.75
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2
## size_pca1   -0.095              
## Trtmnt.rdr2 -0.231  0.006       
## percLoad    -0.954  0.099  0.060
# renaming model to simplify later typing
m2 <- m2.2

##### LRTs for main effects
## Treatment Order
m3 <- update(m2, .~. - Treatment.order)
anova(m2, m3, test = "Chi")
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m3: wbf.aud. ~ size_pca1 + percLoad + (1 | Bee.ID)
## m2: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## m3  5 384.96 395.43 -187.48   374.96                            
## m2  6 378.98 391.54 -183.49   366.98 7.9842      1   0.004719 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT for size
m4 <- update(m2, .~. - size_pca1)
anova(m2, m4)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m4: wbf.aud. ~ Treatment.order + percLoad + (1 | Bee.ID)
## m2: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m4  5 399.76 410.23 -194.88   389.76                             
## m2  6 378.98 391.54 -183.49   366.98 22.783      1  1.813e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# LRT for Treatment (load)
m5 <- update(m2, .~. - percLoad)
anova(m2, m5)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m5: wbf.aud. ~ size_pca1 + Treatment.order + (1 | Bee.ID)
## m2: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## m5  5 383.97 394.45 -186.99   373.97                            
## m2  6 378.98 391.54 -183.49   366.98 6.9967      1   0.008166 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# summarize final model for paper
summary(m2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 360
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.64247 -0.64283  0.06575  0.64332  1.72090 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept)  9.874   3.142   
##  Residual             20.221   4.497   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)      171.9557     3.3390   51.50
## size_pca1         -2.7663     0.4908   -5.64
## Treatment.order2  -3.3999     1.1632   -2.92
## percLoad           5.4191     1.9692    2.75
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2
## size_pca1   -0.095              
## Trtmnt.rdr2 -0.231  0.006       
## percLoad    -0.954  0.099  0.060
# write output
summary(m2)$coefficients  
##                    Estimate Std. Error   t value
## (Intercept)      171.955651  3.3390335 51.498630
## size_pca1         -2.766273  0.4908262 -5.635952
## Treatment.order2  -3.399880  1.1631609 -2.922966
## percLoad           5.419093  1.9691779  2.751957
write.csv( summary(m2)$coefficients, file = "FreqCoefs_percLoad.csv" )


# diagnostics
# qq plot
qqnorm(resid(m2), main = "")
qqline(resid(m2)) # good

# residual plot
plot(fitted(m2), resid(m2), xlab = "fitted", ylab = "residuals")
abline(0,0)

# QQPlot for group-level effects
qqnorm(ranef(m2)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m2)$Bee.ID[[1]]) # looks good

Stroke Amplitude

# make a full model with all two-way interactions
m1 <- lmer(stroke.amplitude ~  (size_pca1 +  Treatment.order + percLoad)^2 + (1|Bee.ID), data = bdta)

summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ (size_pca1 + Treatment.order + percLoad)^2 +  
##     (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 359.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.89733 -0.52961 -0.01694  0.40969  1.99945 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 14.12    3.758   
##  Residual             23.72    4.870   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                            Estimate Std. Error t value
## (Intercept)                 61.9443     6.3305   9.785
## size_pca1                    6.0503     2.3672   2.556
## Treatment.order2             7.4111     9.2790   0.799
## percLoad                    38.4339     3.8380  10.014
## size_pca1:Treatment.order2  -1.1440     0.7947  -1.440
## size_pca1:percLoad          -2.6118     1.3601  -1.920
## Treatment.order2:percLoad   -5.0714     5.7204  -0.887
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 percLd s_1:T. sz_1:L
## size_pca1   -0.201                                   
## Trtmnt.rdr2 -0.825  0.229                            
## percLoad    -0.984  0.176  0.822                     
## sz_pc1:Tr.2  0.259 -0.377 -0.230 -0.258              
## sz_pc1:prcL  0.129 -0.956 -0.172 -0.102  0.217       
## Trtmnt.r2:L  0.814 -0.212 -0.991 -0.826  0.227  0.154
### LRT's for interactions
m2.0 <- update(m1, .~. - Treatment.order:percLoad )

anova(m1, m2.0)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.0: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 | 
## m2.0:     Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## m1: stroke.amplitude ~ (size_pca1 + Treatment.order + percLoad)^2 + 
## m1:     (1 | Bee.ID)
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2.0  8 393.14 409.89 -188.57   377.14                         
## m1    9 394.31 413.16 -188.16   376.31 0.8251      1     0.3637
summary(m2.0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |  
##     Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
##    Data: bdta
## 
## REML criterion at convergence: 365.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.89468 -0.52361 -0.02085  0.43365  2.07383 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 14.75    3.841   
##  Residual             23.18    4.814   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                            Estimate Std. Error t value
## (Intercept)                 66.4391     3.6420  18.242
## size_pca1                    5.5979     2.2913   2.443
## Treatment.order2            -0.7365     1.2581  -0.585
## percLoad                    35.6715     2.1438  16.639
## size_pca1:Treatment.order2  -0.9857     0.7651  -1.288
## size_pca1:percLoad          -2.4201     1.3296  -1.820
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 percLd s_1:T.
## size_pca1   -0.050                            
## Trtmnt.rdr2 -0.227  0.138                     
## percLoad    -0.950  0.002  0.054              
## sz_pc1:Tr.2  0.132 -0.345 -0.035 -0.128       
## sz_pc1:prcL  0.006 -0.955 -0.142  0.046  0.189
m2.1 <- update(m2.0,.~. - size_pca1:Treatment.order)
anova(m2.0, m2.1)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.1: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 | 
## m2.1:     Bee.ID) + size_pca1:percLoad
## m2.0: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 | 
## m2.0:     Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2.1  7 392.99 407.65 -189.49   378.99                         
## m2.0  8 393.14 409.89 -188.57   377.14 1.8488      1     0.1739
summary(m2.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |  
##     Bee.ID) + size_pca1:percLoad
##    Data: bdta
## 
## REML criterion at convergence: 368.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.92050 -0.50058 -0.06054  0.48999  2.08830 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 14.24    3.773   
##  Residual             23.88    4.887   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)         67.1284     3.6558  18.362
## size_pca1            4.5889     2.1787   2.106
## Treatment.order2    -0.7942     1.2763  -0.622
## percLoad            35.2724     2.1552  16.367
## size_pca1:percLoad  -2.1023     1.3241  -1.588
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 percLd
## size_pca1   -0.005                     
## Trtmnt.rdr2 -0.225  0.135              
## percLoad    -0.950 -0.046  0.050       
## sz_pc1:prcL -0.019 -0.966 -0.138  0.072
m2.2 <- update(m2.1, .~. - size_pca1:percLoad)
anova(m2.1, m2.2) ## drop all interactions
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.2: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 | 
## m2.2:     Bee.ID)
## m2.1: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 | 
## m2.1:     Bee.ID) + size_pca1:percLoad
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2.2  6 393.68 406.24 -190.84   381.68                         
## m2.1  7 392.99 407.65 -189.49   378.99 2.6909      1     0.1009
summary(m2.2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |  
##     Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 373.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.87059 -0.50457 -0.03627  0.58077  2.02395 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 14.31    3.782   
##  Residual             24.73    4.973   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)       67.0491     3.7152  18.047
## size_pca1          1.2463     0.5668   2.199
## Treatment.order2  -1.0742     1.2864  -0.835
## percLoad          35.5004     2.1858  16.241
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2
## size_pca1   -0.091              
## Trtmnt.rdr2 -0.230  0.006       
## percLoad    -0.952  0.096  0.060
# renaming model to simplify later typing
m2 <- m2.2

##### LRTs for main effects
## size
m3 <- update(m2, .~. - size_pca1)
anova(m2, m3, test = "Chi") # keep size_pca1
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m3: stroke.amplitude ~ Treatment.order + percLoad + (1 | Bee.ID)
## m2: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 | 
## m2:     Bee.ID)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m3  5 396.50 406.97 -193.25   386.50                           
## m2  6 393.68 406.24 -190.84   381.68 4.8203      1    0.02813 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ Treatment.order + percLoad + (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 379
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.79024 -0.53442  0.00173  0.52524  1.85491 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 17.63    4.199   
##  Residual             24.87    4.987   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)        67.501      3.741   18.04
## Treatment.order2   -1.084      1.290   -0.84
## percLoad           35.221      2.193   16.06
## 
## Correlation of Fixed Effects:
##             (Intr) Trtm.2
## Trtmnt.rdr2 -0.229       
## percLoad    -0.948  0.060
# trt order
m4 <- update(m2, .~. - Treatment.order)
anova(m2, m4) # drop treatment order
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m4: stroke.amplitude ~ size_pca1 + percLoad + (1 | Bee.ID)
## m2: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 | 
## m2:     Bee.ID)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## m4  5 392.42 402.89 -191.21   382.42                         
## m2  6 393.68 406.24 -190.84   381.68 0.7374      1     0.3905
# LRT for Treatment (load)
m5 <- update(m4, .~. - percLoad)
anova(m4, m5)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m5: stroke.amplitude ~ size_pca1 + (1 | Bee.ID)
## m4: stroke.amplitude ~ size_pca1 + percLoad + (1 | Bee.ID)
##    Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
## m5  4 478.07 486.45 -235.04   470.07                             
## m4  5 392.42 402.89 -191.21   382.42 87.657      1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# summarize final model for paper
summary(m4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pca1 + percLoad + (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 376.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.96941 -0.48117 -0.05909  0.48300  1.94459 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 14.53    3.811   
##  Residual             24.42    4.941   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  66.3057     3.5966  18.436
## size_pca1     1.2495     0.5675   2.202
## percLoad     35.6293     2.1693  16.424
## 
## Correlation of Fixed Effects:
##           (Intr) sz_pc1
## size_pca1 -0.091       
## percLoad  -0.965  0.095
# write output
summary(m4)$coefficients  
##              Estimate Std. Error   t value
## (Intercept) 66.305690  3.5965822 18.435750
## size_pca1    1.249463  0.5674547  2.201872
## percLoad    35.629329  2.1693080 16.424283
write.csv( summary(m4)$coefficients, file = "AmpCoefs_percLoad.csv" )


# diagnostics
# qq plot
qqnorm(resid(m4), main = "")
qqline(resid(m4)) # ok

# residual plot
plot(fitted(m4), resid(m3), xlab = "fitted", ylab = "residuals")
abline(0,0)

# QQPlot for group-level effects
qqnorm(ranef(m4)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m4)$Bee.ID[[1]]) # looks good

Wing Velocity

# make a full model with all two-way interactions
m1 <- lmer(wing.velocity ~  (size_pca1 +  Treatment.order + percLoad)^2 + (1|Bee.ID), data = bdta)

summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ (size_pca1 + Treatment.order + percLoad)^2 +  
##     (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: -78.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.00651 -0.41257 -0.06598  0.54758  2.20129 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.003893 0.06239 
##  Residual             0.005935 0.07704 
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                            Estimate Std. Error t value
## (Intercept)                 2.14983    0.10145  21.190
## size_pca1                  -0.12816    0.03757  -3.411
## Treatment.order2           -0.11249    0.14953  -0.752
## percLoad                   -0.41933    0.06150  -6.819
## size_pca1:Treatment.order2  0.01603    0.01258   1.274
## size_pca1:percLoad          0.05043    0.02155   2.340
## Treatment.order2:percLoad   0.05830    0.09222   0.632
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 percLd s_1:T. sz_1:L
## size_pca1   -0.203                                   
## Trtmnt.rdr2 -0.828  0.230                            
## percLoad    -0.984  0.178  0.826                     
## sz_pc1:Tr.2  0.262 -0.377 -0.233 -0.261              
## sz_pc1:prcL  0.131 -0.954 -0.173 -0.103  0.218       
## Trtmnt.r2:L  0.818 -0.214 -0.991 -0.830  0.231  0.155
### LRT's for interactions
m2.0 <- update(m1, .~. - Treatment.order:percLoad )

anova(m1, m2.0)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.0: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 | 
## m2.0:     Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## m1: wing.velocity ~ (size_pca1 + Treatment.order + percLoad)^2 + 
## m1:     (1 | Bee.ID)
##      Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2.0  8 -103.25 -86.492 59.624  -119.25                         
## m1    9 -101.66 -82.811 59.830  -119.66 0.4129      1     0.5205
summary(m2.0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 |  
##     Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
##    Data: bdta
## 
## REML criterion at convergence: -80.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.07164 -0.43946 -0.00926  0.53663  2.21665 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.003959 0.06292 
##  Residual             0.005800 0.07616 
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                            Estimate Std. Error t value
## (Intercept)                 2.09780    0.05777   36.31
## size_pca1                  -0.12306    0.03633   -3.39
## Treatment.order2           -0.01883    0.01990   -0.95
## percLoad                   -0.38734    0.03397  -11.40
## size_pca1:Treatment.order2  0.01421    0.01210    1.17
## size_pca1:percLoad          0.04829    0.02106    2.29
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 percLd s_1:T.
## size_pca1   -0.050                            
## Trtmnt.rdr2 -0.227  0.138                     
## percLoad    -0.949  0.002  0.054              
## sz_pc1:Tr.2  0.132 -0.345 -0.035 -0.129       
## sz_pc1:prcL  0.006 -0.954 -0.142  0.046  0.190
m2.1 <- update(m2.0,.~. - size_pca1:Treatment.order)
anova(m2.0, m2.1)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.1: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 | 
## m2.1:     Bee.ID) + size_pca1:percLoad
## m2.0: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 | 
## m2.0:     Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
##      Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m2.1  7 -103.70 -89.044 58.852  -117.70                         
## m2.0  8 -103.25 -86.492 59.624  -119.25 1.5428      1     0.2142
summary(m2.1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 |  
##     Bee.ID) + size_pca1:percLoad
##    Data: bdta
## 
## REML criterion at convergence: -86.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.08227 -0.51145 -0.04384  0.65437  2.23563 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.003885 0.06233 
##  Residual             0.005902 0.07682 
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)         2.08849    0.05768   36.21
## size_pca1          -0.10840    0.03436   -3.16
## Treatment.order2   -0.01800    0.02006   -0.90
## percLoad           -0.38198    0.03395  -11.25
## size_pca1:percLoad  0.04363    0.02084    2.09
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 Trtm.2 percLd
## size_pca1   -0.005                     
## Trtmnt.rdr2 -0.224  0.135              
## percLoad    -0.948 -0.046  0.050       
## sz_pc1:prcL -0.019 -0.965 -0.138  0.072
m2.2 <- update(m2.1, .~. - size_pca1:percLoad)
anova(m2.1, m2.2) ## don't drop size:percLoad interaction
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m2.2: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 | 
## m2.2:     Bee.ID)
## m2.1: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 | 
## m2.1:     Bee.ID) + size_pca1:percLoad
##      Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)  
## m2.2  6 -101.16 -88.596 56.581  -113.16                           
## m2.1  7 -103.70 -89.044 58.852  -117.70 4.5421      1    0.03307 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# renaming model to simplify later typing
m2 <- m2.1

##### LRTs for main effects
# trt order
m3 <- update(m2, .~. - Treatment.order)
anova(m2, m3, test = "Chi") # drop trt order
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## m3: wing.velocity ~ size_pca1 + percLoad + (1 | Bee.ID) + size_pca1:percLoad
## m2: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 | 
## m2:     Bee.ID) + size_pca1:percLoad
##    Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## m3  6 -104.83 -92.261 58.413  -116.83                         
## m2  7 -103.70 -89.044 58.852  -117.70 0.8775      1     0.3489
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## wing.velocity ~ size_pca1 + percLoad + (1 | Bee.ID) + size_pca1:percLoad
##    Data: bdta
## 
## REML criterion at convergence: -91.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.18903 -0.53485 -0.00798  0.63629  2.12781 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.003898 0.06243 
##  Residual             0.005865 0.07658 
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)         2.07698    0.05606   37.05
## size_pca1          -0.10424    0.03395   -3.07
## percLoad           -0.38053    0.03381  -11.25
## size_pca1:percLoad  0.04104    0.02058    1.99
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 percLd
## size_pca1    0.026              
## percLoad    -0.963 -0.053       
## sz_pc1:prcL -0.052 -0.964  0.080
# summarize final model for paper
summary(m3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## wing.velocity ~ size_pca1 + percLoad + (1 | Bee.ID) + size_pca1:percLoad
##    Data: bdta
## 
## REML criterion at convergence: -91.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.18903 -0.53485 -0.00798  0.63629  2.12781 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.003898 0.06243 
##  Residual             0.005865 0.07658 
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                    Estimate Std. Error t value
## (Intercept)         2.07698    0.05606   37.05
## size_pca1          -0.10424    0.03395   -3.07
## percLoad           -0.38053    0.03381  -11.25
## size_pca1:percLoad  0.04104    0.02058    1.99
## 
## Correlation of Fixed Effects:
##             (Intr) sz_pc1 percLd
## size_pca1    0.026              
## percLoad    -0.963 -0.053       
## sz_pc1:prcL -0.052 -0.964  0.080
# write output
summary(m3)$coefficients  
##                       Estimate Std. Error    t value
## (Intercept)         2.07697845 0.05606361  37.046819
## size_pca1          -0.10424096 0.03394778  -3.070627
## percLoad           -0.38052624 0.03381059 -11.254646
## size_pca1:percLoad  0.04104281 0.02058223   1.994090
write.csv( summary(m3)$coefficients, file = "VelocCoefs_percLoad.csv" )


# diagnostics
# qq plot
qqnorm(resid(m3), main = "")
qqline(resid(m3)) # good

# residual plot
plot(fitted(m3), resid(m3), xlab = "fitted", ylab = "residuals")
abline(0,0)

# QQPlot for group-level effects
qqnorm(ranef(m3)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m3)$Bee.ID[[1]]) # looks good

Part II: Which flight measurements are most closely associated with

respiratory rate in bees – wingbeat frequency, stroke amplitude, or

wing velocity?

Multilevel model approach to see the relatedness of the response variables

  • Wingbeat frequency
  • Stroke amplitude
  • Wing velocity
  • Resp. Rate
# predicting avg resp rate, including stroke amplitude, wing velocity and wingbeat frequency
# looks like we have to remove wing velocity
vif.mer(lmer(av.resp..CO2.mL.hr. ~ stroke.amplitude + wing.velocity + wbf.aud. + percLoad + Treatment.order + size_pca1 + (1|Bee.ID), data = bdta))
## stroke.amplitude    wing.velocity         wbf.aud.         percLoad 
##       190.642034       133.592345        11.644262        10.853222 
## Treatment.order2        size_pca1 
##         1.262757         1.364694
vif.mer(lmer(av.resp..CO2.mL.hr. ~ stroke.amplitude + wbf.aud. + percLoad + Treatment.order + size_pca1 + (1|Bee.ID), data = bdta)) # still worryingly high for percLoad and stroke amplitude
## stroke.amplitude         wbf.aud.         percLoad Treatment.order2 
##         7.741481         1.830209         8.182520         1.255652 
##        size_pca1 
##         1.355159
mm1 <- lmer(av.resp..CO2.mL.hr. ~ (stroke.amplitude + wbf.aud. + percLoad + Treatment.order + size_pca1)^2 + (1|Bee.ID), data = bdta)
## Warning: Some predictor variables are on very different scales: consider
## rescaling
summary(mm1)  # warning says to rescale
## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ (stroke.amplitude + wbf.aud. + percLoad +  
##     Treatment.order + size_pca1)^2 + (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 211.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4812 -0.3747 -0.1374  0.2959  2.1473 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.9021   0.9498  
##  Residual             0.4039   0.6356  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                     Estimate Std. Error t value
## (Intercept)                       -59.720814  44.421604  -1.344
## stroke.amplitude                    0.359487   0.608891   0.590
## wbf.aud.                            0.278405   0.254698   1.093
## percLoad                           11.733723  29.320882   0.400
## Treatment.order2                    9.248884   8.744342   1.058
## size_pca1                           3.953128   2.787254   1.418
## stroke.amplitude:wbf.aud.          -0.001567   0.003489  -0.449
## stroke.amplitude:percLoad          -0.089640   0.050249  -1.784
## stroke.amplitude:Treatment.order2   0.019955   0.036748   0.543
## stroke.amplitude:size_pca1          0.010180   0.016490   0.617
## wbf.aud.:percLoad                   0.033891   0.149234   0.227
## wbf.aud.:Treatment.order2          -0.047825   0.043855  -1.091
## wbf.aud.:size_pca1                 -0.017492   0.013469  -1.299
## percLoad:Treatment.order2          -1.939776   1.533297  -1.265
## percLoad:size_pca1                 -0.158919   0.670849  -0.237
## Treatment.order2:size_pca1         -0.455025   0.180000  -2.528
## 
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
# rescale numeric predictors
bdta_scaled = scale(bdta[, c("av.resp..CO2.mL.hr.", 
                             "stroke.amplitude", 
                             "wbf.aud.",
                             "percLoad",
                             "size_pca1")], 
                    center = TRUE, 
                    scale = TRUE)
colnames(bdta_scaled) <- paste0(colnames(bdta_scaled), "_scaled")
bdta = cbind(bdta, bdta_scaled)

mm1 <- lmer(av.resp..CO2.mL.hr. ~ (stroke.amplitude_scaled + wbf.aud._scaled + percLoad + Treatment.order + size_pca1_scaled)^2 + (1|Bee.ID), data = bdta)
summary(mm1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## av.resp..CO2.mL.hr. ~ (stroke.amplitude_scaled + wbf.aud._scaled +  
##     percLoad + Treatment.order + size_pca1_scaled)^2 + (1 | Bee.ID)
##    Data: bdta
## 
## REML criterion at convergence: 160.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4812 -0.3747 -0.1374  0.2959  2.1473 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.9021   0.9498  
##  Residual             0.4039   0.6356  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                          Estimate Std. Error t value
## (Intercept)                               -0.1472     2.1940  -0.067
## stroke.amplitude_scaled                    0.9720     0.9560   1.017
## wbf.aud._scaled                            0.6494     1.7496   0.371
## percLoad                                   6.7446     1.4013   4.813
## Treatment.order2                           3.1524     2.4584   1.282
## size_pca1_scaled                           3.5029     1.7806   1.967
## stroke.amplitude_scaled:wbf.aud._scaled   -0.1466     0.3266  -0.449
## stroke.amplitude_scaled:percLoad          -1.1008     0.6171  -1.784
## stroke.amplitude_scaled:Treatment.order2   0.2451     0.4513   0.543
## stroke.amplitude_scaled:size_pca1_scaled   0.2107     0.3413   0.617
## wbf.aud._scaled:percLoad                   0.2583     1.1374   0.227
## wbf.aud._scaled:Treatment.order2          -0.3645     0.3342  -1.091
## wbf.aud._scaled:size_pca1_scaled          -0.2247     0.1730  -1.299
## percLoad:Treatment.order2                 -1.9398     1.5333  -1.265
## percLoad:size_pca1_scaled                 -0.2678     1.1305  -0.237
## Treatment.order2:size_pca1_scaled         -0.7668     0.3033  -2.528
## 
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
# LRT's
mm2 <- update(mm1, .~. - wbf.aud._scaled:size_pca1_scaled)
anova(mm1, mm2)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm2: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm2:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm2:     stroke.amplitude_scaled:wbf.aud._scaled + stroke.amplitude_scaled:percLoad + 
## mm2:     stroke.amplitude_scaled:Treatment.order + stroke.amplitude_scaled:size_pca1_scaled + 
## mm2:     wbf.aud._scaled:percLoad + wbf.aud._scaled:Treatment.order + 
## mm2:     percLoad:Treatment.order + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
## mm1: av.resp..CO2.mL.hr. ~ (stroke.amplitude_scaled + wbf.aud._scaled + 
## mm1:     percLoad + Treatment.order + size_pca1_scaled)^2 + (1 | Bee.ID)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mm2 17 183.70 219.31 -74.852   149.70                         
## mm1 18 183.48 221.18 -73.741   147.48 2.2212      1     0.1361
summary(mm2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +  
##     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +  
##     stroke.amplitude_scaled:wbf.aud._scaled + stroke.amplitude_scaled:percLoad +  
##     stroke.amplitude_scaled:Treatment.order + stroke.amplitude_scaled:size_pca1_scaled +  
##     wbf.aud._scaled:percLoad + wbf.aud._scaled:Treatment.order +  
##     percLoad:Treatment.order + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 160.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.49191 -0.41271 -0.08384  0.31098  2.15506 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.9109   0.9544  
##  Residual             0.4120   0.6419  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                          Estimate Std. Error t value
## (Intercept)                                0.3898     2.1726   0.179
## stroke.amplitude_scaled                    0.9845     0.9643   1.021
## wbf.aud._scaled                            0.1441     1.7214   0.084
## percLoad                                   6.4778     1.3979   4.634
## Treatment.order2                           2.8888     2.4719   1.169
## size_pca1_scaled                           3.8549     1.7749   2.172
## stroke.amplitude_scaled:wbf.aud._scaled   -0.2753     0.3142  -0.876
## stroke.amplitude_scaled:percLoad          -1.0519     0.6214  -1.693
## stroke.amplitude_scaled:Treatment.order2   0.1989     0.4541   0.438
## stroke.amplitude_scaled:size_pca1_scaled   0.2653     0.3414   0.777
## wbf.aud._scaled:percLoad                   0.5750     1.1210   0.513
## wbf.aud._scaled:Treatment.order2          -0.3539     0.3373  -1.049
## percLoad:Treatment.order2                 -1.7899     1.5426  -1.160
## percLoad:size_pca1_scaled                 -0.5482     1.1190  -0.490
## Treatment.order2:size_pca1_scaled         -0.6698     0.2968  -2.257
## 
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
mm3 <- update(mm2, .~. -Treatment.order:size_pca1_scaled) # keep interaction
anova(mm3, mm2)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm3: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm3:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm3:     stroke.amplitude_scaled:wbf.aud._scaled + stroke.amplitude_scaled:percLoad + 
## mm3:     stroke.amplitude_scaled:Treatment.order + stroke.amplitude_scaled:size_pca1_scaled + 
## mm3:     wbf.aud._scaled:percLoad + wbf.aud._scaled:Treatment.order + 
## mm3:     percLoad:Treatment.order + percLoad:size_pca1_scaled
## mm2: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm2:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm2:     stroke.amplitude_scaled:wbf.aud._scaled + stroke.amplitude_scaled:percLoad + 
## mm2:     stroke.amplitude_scaled:Treatment.order + stroke.amplitude_scaled:size_pca1_scaled + 
## mm2:     wbf.aud._scaled:percLoad + wbf.aud._scaled:Treatment.order + 
## mm2:     percLoad:Treatment.order + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
## mm3 16 188.86 222.37 -78.431   156.86                            
## mm2 17 183.70 219.31 -74.852   149.70 7.1571      1   0.007467 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mm3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +  
##     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +  
##     stroke.amplitude_scaled:wbf.aud._scaled + stroke.amplitude_scaled:percLoad +  
##     stroke.amplitude_scaled:Treatment.order + stroke.amplitude_scaled:size_pca1_scaled +  
##     wbf.aud._scaled:percLoad + wbf.aud._scaled:Treatment.order +  
##     percLoad:Treatment.order + percLoad:size_pca1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 164.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3975 -0.5111 -0.0859  0.4055  2.1634 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.8748   0.9353  
##  Residual             0.4932   0.7023  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                          Estimate Std. Error t value
## (Intercept)                                1.4031     2.2484   0.624
## stroke.amplitude_scaled                    0.4383     1.0037   0.437
## wbf.aud._scaled                            0.1332     1.8467   0.072
## percLoad                                   5.8215     1.4465   4.025
## Treatment.order2                           3.0253     2.6625   1.136
## size_pca1_scaled                           3.7037     1.8855   1.964
## stroke.amplitude_scaled:wbf.aud._scaled   -0.3434     0.3355  -1.024
## stroke.amplitude_scaled:percLoad          -0.5762     0.6295  -0.915
## stroke.amplitude_scaled:Treatment.order2   0.1250     0.4886   0.256
## stroke.amplitude_scaled:size_pca1_scaled   0.2303     0.3602   0.639
## wbf.aud._scaled:percLoad                   0.3812     1.1934   0.319
## wbf.aud._scaled:Treatment.order2           0.1745     0.2655   0.657
## percLoad:Treatment.order2                 -1.8886     1.6610  -1.137
## percLoad:size_pca1_scaled                 -0.7152     1.1828  -0.605
## 
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
mm4 <- update(mm2, .~. - wbf.aud._scaled:Treatment.order)
anova(mm2, mm4)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm4: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm4:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm4:     stroke.amplitude_scaled:wbf.aud._scaled + stroke.amplitude_scaled:percLoad + 
## mm4:     stroke.amplitude_scaled:Treatment.order + stroke.amplitude_scaled:size_pca1_scaled + 
## mm4:     wbf.aud._scaled:percLoad + percLoad:Treatment.order + percLoad:size_pca1_scaled + 
## mm4:     Treatment.order:size_pca1_scaled
## mm2: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm2:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm2:     stroke.amplitude_scaled:wbf.aud._scaled + stroke.amplitude_scaled:percLoad + 
## mm2:     stroke.amplitude_scaled:Treatment.order + stroke.amplitude_scaled:size_pca1_scaled + 
## mm2:     wbf.aud._scaled:percLoad + wbf.aud._scaled:Treatment.order + 
## mm2:     percLoad:Treatment.order + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mm4 16 183.54 217.05 -75.769   151.54                         
## mm2 17 183.70 219.31 -74.852   149.70 1.8345      1     0.1756
summary(mm4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +  
##     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +  
##     stroke.amplitude_scaled:wbf.aud._scaled + stroke.amplitude_scaled:percLoad +  
##     stroke.amplitude_scaled:Treatment.order + stroke.amplitude_scaled:size_pca1_scaled +  
##     wbf.aud._scaled:percLoad + percLoad:Treatment.order + percLoad:size_pca1_scaled +  
##     Treatment.order:size_pca1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 161.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.54564 -0.45657 -0.08531  0.37288  2.24119 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.8871   0.9418  
##  Residual             0.4227   0.6502  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                          Estimate Std. Error t value
## (Intercept)                               0.09727    2.17454   0.045
## stroke.amplitude_scaled                   0.67101    0.92807   0.723
## wbf.aud._scaled                           0.01824    1.73163   0.011
## percLoad                                  6.66930    1.39819   4.770
## Treatment.order2                          3.77686    2.36346   1.598
## size_pca1_scaled                          3.69859    1.77947   2.078
## stroke.amplitude_scaled:wbf.aud._scaled  -0.29890    0.31607  -0.946
## stroke.amplitude_scaled:percLoad         -0.88041    0.60607  -1.453
## stroke.amplitude_scaled:Treatment.order2  0.29973    0.44709   0.670
## stroke.amplitude_scaled:size_pca1_scaled  0.19508    0.33683   0.579
## wbf.aud._scaled:percLoad                  0.52625    1.12791   0.467
## percLoad:Treatment.order2                -2.35292    1.47151  -1.599
## percLoad:size_pca1_scaled                -0.53345    1.12492  -0.474
## Treatment.order2:size_pca1_scaled        -0.45825    0.21875  -2.095
## 
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
mm5 <- update(mm4, .~. - stroke.amplitude_scaled:wbf.aud._scaled)
anova(mm4, mm5)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm5: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm5:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm5:     stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:Treatment.order + 
## mm5:     stroke.amplitude_scaled:size_pca1_scaled + wbf.aud._scaled:percLoad + 
## mm5:     percLoad:Treatment.order + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
## mm4: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm4:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm4:     stroke.amplitude_scaled:wbf.aud._scaled + stroke.amplitude_scaled:percLoad + 
## mm4:     stroke.amplitude_scaled:Treatment.order + stroke.amplitude_scaled:size_pca1_scaled + 
## mm4:     wbf.aud._scaled:percLoad + percLoad:Treatment.order + percLoad:size_pca1_scaled + 
## mm4:     Treatment.order:size_pca1_scaled
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mm5 15 182.56 213.98 -76.281   152.56                         
## mm4 16 183.54 217.05 -75.769   151.54 1.0242      1     0.3115
summary(mm5)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +  
##     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +  
##     stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:Treatment.order +  
##     stroke.amplitude_scaled:size_pca1_scaled + wbf.aud._scaled:percLoad +  
##     percLoad:Treatment.order + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 161.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6390 -0.4133 -0.1268  0.3654  2.3371 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.9212   0.9598  
##  Residual             0.4080   0.6388  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                          Estimate Std. Error t value
## (Intercept)                                0.1984     2.1527   0.092
## stroke.amplitude_scaled                    0.9261     0.8850   1.047
## wbf.aud._scaled                            1.4238     0.8461   1.683
## percLoad                                   6.6563     1.3848   4.807
## Treatment.order2                           4.0182     2.3080   1.741
## size_pca1_scaled                           4.4330     1.5838   2.799
## stroke.amplitude_scaled:percLoad          -1.0551     0.5757  -1.833
## stroke.amplitude_scaled:Treatment.order2   0.3835     0.4335   0.885
## stroke.amplitude_scaled:size_pca1_scaled   0.3313     0.3041   1.089
## wbf.aud._scaled:percLoad                  -0.3826     0.5739  -0.667
## percLoad:Treatment.order2                 -2.5079     1.4368  -1.745
## percLoad:size_pca1_scaled                 -1.0040     1.0016  -1.002
## Treatment.order2:size_pca1_scaled         -0.4638     0.2153  -2.154
## 
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it
mm6 <- update(mm5, .~. - stroke.amplitude_scaled:Treatment.order)
anova(mm5, mm6)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm6: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm6:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm6:     stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled + 
## mm6:     wbf.aud._scaled:percLoad + percLoad:Treatment.order + percLoad:size_pca1_scaled + 
## mm6:     Treatment.order:size_pca1_scaled
## mm5: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm5:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm5:     stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:Treatment.order + 
## mm5:     stroke.amplitude_scaled:size_pca1_scaled + wbf.aud._scaled:percLoad + 
## mm5:     percLoad:Treatment.order + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mm6 14 181.81 211.13 -76.903   153.81                         
## mm5 15 182.56 213.98 -76.281   152.56 1.2439      1     0.2647
summary(mm6)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +  
##     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +  
##     stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled +  
##     wbf.aud._scaled:percLoad + percLoad:Treatment.order + percLoad:size_pca1_scaled +  
##     Treatment.order:size_pca1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 162.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.67475 -0.46717 -0.06738  0.39201  2.38883 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.8984   0.9479  
##  Residual             0.4128   0.6425  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                          Estimate Std. Error t value
## (Intercept)                                0.9357     1.9719   0.474
## stroke.amplitude_scaled                    1.2688     0.7898   1.607
## wbf.aud._scaled                            1.5030     0.8458   1.777
## percLoad                                   6.2134     1.2846   4.837
## Treatment.order2                           2.7598     1.7714   1.558
## size_pca1_scaled                           4.3302     1.5810   2.739
## stroke.amplitude_scaled:percLoad          -1.1420     0.5668  -2.015
## stroke.amplitude_scaled:size_pca1_scaled   0.2958     0.3017   0.981
## wbf.aud._scaled:percLoad                  -0.4589     0.5695  -0.806
## percLoad:Treatment.order2                 -1.7307     1.1093  -1.560
## percLoad:size_pca1_scaled                 -0.9794     1.0021  -0.977
## Treatment.order2:size_pca1_scaled         -0.4111     0.2072  -1.984
## 
## Correlation of Fixed Effects:
##             (Intr) strk._ wbf.._ percLd Trtm.2 sz_p1_ st._:L s._:_1 w.._:L
## strk.mpltd_ -0.223                                                        
## wbf.d._scld -0.060 -0.037                                                 
## percLoad    -0.991  0.283  0.065                                          
## Trtmnt.rdr2 -0.545  0.057  0.375  0.530                                   
## sz_pc1_scld -0.191  0.198  0.441  0.174  0.153                            
## strk.mpl_:L  0.485 -0.939  0.062 -0.542 -0.082 -0.207                     
## strk.m_:_1_  0.005  0.214  0.080 -0.022 -0.105  0.841 -0.174              
## wbf.d._sc:L  0.161 -0.087 -0.967 -0.178 -0.396 -0.427  0.075 -0.072       
## prcLd:Trt.2  0.562 -0.070 -0.370 -0.554 -0.994 -0.132  0.101  0.125  0.405
## prcLd:sz_1_  0.170 -0.236 -0.461 -0.158 -0.138 -0.985  0.231 -0.852  0.469
## Trtmn.2:_1_  0.378 -0.161  0.198 -0.381 -0.218  0.094  0.258  0.240 -0.190
##             pL:T.2 pL:_1_
## strk.mpltd_              
## wbf.d._scld              
## percLoad                 
## Trtmnt.rdr2              
## sz_pc1_scld              
## strk.mpl_:L              
## strk.m_:_1_              
## wbf.d._sc:L              
## prcLd:Trt.2              
## prcLd:sz_1_  0.122       
## Trtmn.2:_1_  0.219 -0.171
mm7 <- update(mm6, .~. - percLoad:Treatment.order)
anova(mm6, mm7)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm7: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm7:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm7:     stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled + 
## mm7:     wbf.aud._scaled:percLoad + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
## mm6: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm6:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm6:     stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled + 
## mm6:     wbf.aud._scaled:percLoad + percLoad:Treatment.order + percLoad:size_pca1_scaled + 
## mm6:     Treatment.order:size_pca1_scaled
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mm7 13 181.80 209.02 -77.899   155.80                         
## mm6 14 181.81 211.13 -76.903   153.81 1.9911      1     0.1582
summary(mm7)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +  
##     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +  
##     stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled +  
##     wbf.aud._scaled:percLoad + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 167
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.71405 -0.37038 -0.02943  0.33275  2.39389 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.0729   1.0358  
##  Residual             0.3714   0.6094  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                           Estimate Std. Error t value
## (Intercept)                               2.634514   1.608837   1.638
## stroke.amplitude_scaled                   1.315040   0.769176   1.710
## wbf.aud._scaled                           1.049892   0.754406   1.392
## percLoad                                  5.141776   1.055823   4.870
## Treatment.order2                         -0.002183   0.188595  -0.012
## size_pca1_scaled                          4.065746   1.542299   2.636
## stroke.amplitude_scaled:percLoad         -1.137541   0.551286  -2.063
## stroke.amplitude_scaled:size_pca1_scaled  0.361927   0.296016   1.223
## wbf.aud._scaled:percLoad                 -0.150894   0.503708  -0.300
## percLoad:size_pca1_scaled                -0.850048   0.981109  -0.866
## Treatment.order2:size_pca1_scaled        -0.336078   0.193590  -1.736
## 
## Correlation of Fixed Effects:
##             (Intr) strk._ wbf.._ percLd Trtm.2 sz_p1_ st._:L s._:_1 w.._:L
## strk.mpltd_ -0.226                                                        
## wbf.d._scld  0.179 -0.049                                                 
## percLoad    -0.987  0.297 -0.168                                          
## Trtmnt.rdr2  0.165 -0.132  0.069 -0.237                                   
## sz_pc1_scld -0.126  0.212  0.432  0.107  0.194                            
## strk.mpl_:L  0.528 -0.938  0.089 -0.593  0.188 -0.207                     
## strk.m_:_1_ -0.053  0.242  0.152  0.032  0.169  0.877 -0.192              
## wbf.d._sc:L -0.067 -0.094 -0.959  0.039  0.072 -0.421  0.068 -0.154       
## prcLd:sz_1_  0.109 -0.253 -0.458 -0.097 -0.151 -0.984  0.235 -0.887  0.472
## Trtmn.2:_1_  0.324 -0.134  0.311 -0.327 -0.006  0.145  0.235  0.239 -0.314
##             pL:_1_
## strk.mpltd_       
## wbf.d._scld       
## percLoad          
## Trtmnt.rdr2       
## sz_pc1_scld       
## strk.mpl_:L       
## strk.m_:_1_       
## wbf.d._sc:L       
## prcLd:sz_1_       
## Trtmn.2:_1_ -0.221
mm8 <- update(mm7, .~. - wbf.aud._scaled:percLoad)
anova(mm8, mm7)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm8: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm8:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm8:     stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled + 
## mm8:     percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
## mm7: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm7:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm7:     stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled + 
## mm7:     wbf.aud._scaled:percLoad + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mm8 12 179.99 205.12 -77.996   155.99                         
## mm7 13 181.80 209.02 -77.899   155.80 0.1943      1     0.6594
summary(mm8)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +  
##     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +  
##     stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled +  
##     percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 167.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.73514 -0.34367 -0.05258  0.32887  2.40909 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.0586   1.0289  
##  Residual             0.3631   0.6026  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                           Estimate Std. Error t value
## (Intercept)                               2.600949   1.589185   1.637
## stroke.amplitude_scaled                   1.297131   0.757782   1.712
## wbf.aud._scaled                           0.831911   0.211132   3.940
## percLoad                                  5.155421   1.044492   4.936
## Treatment.order2                          0.001467   0.186058   0.008
## size_pca1_scaled                          3.871007   1.384576   2.796
## stroke.amplitude_scaled:percLoad         -1.128710   0.544304  -2.074
## stroke.amplitude_scaled:size_pca1_scaled  0.348367   0.289558   1.203
## percLoad:size_pca1_scaled                -0.711749   0.856171  -0.831
## Treatment.order2:size_pca1_scaled        -0.354355   0.181766  -1.950
## 
## Correlation of Fixed Effects:
##             (Intr) strk._ wbf.._ percLd Trtm.2 sz_p1_ st._:L s._:_1 pL:_1_
## strk.mpltd_ -0.234                                                        
## wbf.d._scld  0.407 -0.497                                                 
## percLoad    -0.987  0.302 -0.462                                          
## Trtmnt.rdr2  0.172 -0.126  0.492 -0.241                                   
## sz_pc1_scld -0.169  0.191  0.107  0.136  0.248                            
## strk.mpl_:L  0.535 -0.938  0.549 -0.598  0.184 -0.196                     
## strk.m_:_1_ -0.063  0.231  0.015  0.037  0.183  0.906 -0.184              
## prcLd:sz_1_  0.160 -0.238 -0.020 -0.131 -0.211 -0.982  0.231 -0.935       
## Trtmn.2:_1_  0.320 -0.173  0.034 -0.332  0.018  0.015  0.270  0.204 -0.087
mm9 <- update(mm8, .~. - percLoad:size_pca1_scaled)
anova(mm9, mm8)
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm9: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm9:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm9:     stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled + 
## mm9:     Treatment.order:size_pca1_scaled
## mm8: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm8:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm8:     stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled + 
## mm8:     percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mm9 11 178.87 201.91 -78.436   156.87                         
## mm8 12 179.99 205.12 -77.996   155.99 0.8806      1      0.348
summary(mm9)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +  
##     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +  
##     stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled +  
##     Treatment.order:size_pca1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 169.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.87306 -0.37383 -0.07975  0.32076  2.55486 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.0512   1.025   
##  Residual             0.3612   0.601   
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                          Estimate Std. Error t value
## (Intercept)                               2.81219    1.56426   1.798
## stroke.amplitude_scaled                   1.14638    0.73391   1.562
## wbf.aud._scaled                           0.82870    0.21048   3.937
## percLoad                                  5.04157    1.03253   4.883
## Treatment.order2                         -0.03107    0.18138  -0.171
## size_pca1_scaled                          2.74118    0.26284  10.429
## stroke.amplitude_scaled:percLoad         -1.02381    0.52811  -1.939
## stroke.amplitude_scaled:size_pca1_scaled  0.12336    0.10254   1.203
## Treatment.order2:size_pca1_scaled        -0.36752    0.18058  -2.035
## 
## Correlation of Fixed Effects:
##             (Intr) strk._ wbf.._ percLd Trtm.2 sz_p1_ st._:L s._:_1
## strk.mpltd_ -0.204                                                 
## wbf.d._scld  0.415 -0.516                                          
## percLoad    -0.987  0.282 -0.468                                   
## Trtmnt.rdr2  0.213 -0.186  0.499 -0.277                            
## sz_pc1_scld -0.068 -0.234  0.463  0.039  0.222                     
## strk.mpl_:L  0.518 -0.934  0.569 -0.589  0.245  0.163              
## strk.m_:_1_  0.245  0.025 -0.008 -0.241 -0.042 -0.170  0.091       
## Trtmn.2:_1_  0.340 -0.200  0.033 -0.348  0.000 -0.374  0.300  0.345
mm10 <- update(mm9, .~. - stroke.amplitude_scaled:percLoad)
anova(mm10, mm9) # keep interaction
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm10: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm10:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm10:     stroke.amplitude_scaled:size_pca1_scaled + Treatment.order:size_pca1_scaled
## mm9: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm9:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm9:     stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled + 
## mm9:     Treatment.order:size_pca1_scaled
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## mm10 10 181.17 202.12 -80.587   161.17                           
## mm9  11 178.87 201.91 -78.436   156.87 4.3005      1     0.0381 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mm10)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +  
##     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +  
##     stroke.amplitude_scaled:size_pca1_scaled + Treatment.order:size_pca1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 173.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.92324 -0.37734 -0.04338  0.31256  2.70415 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 0.9448   0.9720  
##  Residual             0.4257   0.6524  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                          Estimate Std. Error t value
## (Intercept)                               4.27791    1.40239   3.050
## stroke.amplitude_scaled                  -0.20299    0.27398  -0.741
## wbf.aud._scaled                           1.06817    0.18262   5.849
## percLoad                                  3.92936    0.87680   4.481
## Treatment.order2                          0.05687    0.18987   0.300
## size_pca1_scaled                          2.83452    0.26014  10.896
## stroke.amplitude_scaled:size_pca1_scaled  0.14257    0.10966   1.300
## Treatment.order2:size_pca1_scaled        -0.26620    0.18609  -1.431
## 
## Correlation of Fixed Effects:
##             (Intr) strk._ wbf.._ percLd Trtm.2 sz_p1_ s._:_1
## strk.mpltd_  0.918                                          
## wbf.d._scld  0.183  0.066                                   
## percLoad    -0.988 -0.927 -0.213                            
## Trtmnt.rdr2  0.102  0.129  0.442 -0.170                     
## sz_pc1_scld -0.182 -0.233  0.476  0.169  0.195              
## strk.m_:_1_  0.214  0.291 -0.061 -0.214 -0.064 -0.189       
## Trtmn.2:_1_  0.215  0.222 -0.167 -0.211 -0.076 -0.474  0.326
mm11 <- update(mm9, .~. - stroke.amplitude_scaled:size_pca1_scaled)
anova(mm9, mm11) # drop interaction
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm11: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm11:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm11:     stroke.amplitude_scaled:percLoad + Treatment.order:size_pca1_scaled
## mm9: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm9:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm9:     stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled + 
## mm9:     Treatment.order:size_pca1_scaled
##      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mm11 10 178.60 199.54 -79.300   158.60                         
## mm9  11 178.87 201.91 -78.436   156.87 1.7279      1     0.1887
summary(mm11)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +  
##     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +  
##     stroke.amplitude_scaled:percLoad + Treatment.order:size_pca1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 168.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.75223 -0.41953 -0.02529  0.44228  2.43853 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.0672   1.0331  
##  Residual             0.3623   0.6019  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                   Estimate Std. Error t value
## (Intercept)                        2.34652    1.52140   1.542
## stroke.amplitude_scaled            1.12885    0.73555   1.535
## wbf.aud._scaled                    0.82929    0.21120   3.927
## percLoad                           5.34378    1.00541   5.315
## Treatment.order2                  -0.02245    0.18159  -0.124
## size_pca1_scaled                   2.79423    0.26035  10.732
## stroke.amplitude_scaled:percLoad  -1.08484    0.52736  -2.057
## Treatment.order2:size_pca1_scaled -0.44257    0.16981  -2.606
## 
## Correlation of Fixed Effects:
##             (Intr) strk._ wbf.._ percLd Trtm.2 sz_p1_ st._:L
## strk.mpltd_ -0.217                                          
## wbf.d._scld  0.431 -0.517                                   
## percLoad    -0.986  0.297 -0.485                            
## Trtmnt.rdr2  0.231 -0.186  0.500 -0.296                     
## sz_pc1_scld -0.027 -0.233  0.468 -0.002  0.218              
## strk.mpl_:L  0.514 -0.941  0.573 -0.587  0.251  0.182       
## Trtmn.2:_1_  0.281 -0.222  0.038 -0.291  0.016 -0.339  0.287
## main effects
mm14 <- update(mm11, .~. - wbf.aud._scaled)
anova(mm11, mm14) # wbf is significant
## refitting model(s) with ML (instead of REML)
## Data: bdta
## Models:
## mm14: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + percLoad + Treatment.order + 
## mm14:     size_pca1_scaled + (1 | Bee.ID) + stroke.amplitude_scaled:percLoad + 
## mm14:     Treatment.order:size_pca1_scaled
## mm11: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled + 
## mm11:     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) + 
## mm11:     stroke.amplitude_scaled:percLoad + Treatment.order:size_pca1_scaled
##      Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)    
## mm14  9 191.71 210.56 -86.855   173.71                            
## mm11 10 178.60 199.54 -79.300   158.60 15.11      1  0.0001014 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mm11)
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +  
##     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +  
##     stroke.amplitude_scaled:percLoad + Treatment.order:size_pca1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 168.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.75223 -0.41953 -0.02529  0.44228  2.43853 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.0672   1.0331  
##  Residual             0.3623   0.6019  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                   Estimate Std. Error t value
## (Intercept)                        2.34652    1.52140   1.542
## stroke.amplitude_scaled            1.12885    0.73555   1.535
## wbf.aud._scaled                    0.82929    0.21120   3.927
## percLoad                           5.34378    1.00541   5.315
## Treatment.order2                  -0.02245    0.18159  -0.124
## size_pca1_scaled                   2.79423    0.26035  10.732
## stroke.amplitude_scaled:percLoad  -1.08484    0.52736  -2.057
## Treatment.order2:size_pca1_scaled -0.44257    0.16981  -2.606
## 
## Correlation of Fixed Effects:
##             (Intr) strk._ wbf.._ percLd Trtm.2 sz_p1_ st._:L
## strk.mpltd_ -0.217                                          
## wbf.d._scld  0.431 -0.517                                   
## percLoad    -0.986  0.297 -0.485                            
## Trtmnt.rdr2  0.231 -0.186  0.500 -0.296                     
## sz_pc1_scld -0.027 -0.233  0.468 -0.002  0.218              
## strk.mpl_:L  0.514 -0.941  0.573 -0.587  0.251  0.182       
## Trtmn.2:_1_  0.281 -0.222  0.038 -0.291  0.016 -0.339  0.287
summary(mm11)  # final model for paper
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +  
##     percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +  
##     stroke.amplitude_scaled:percLoad + Treatment.order:size_pca1_scaled
##    Data: bdta
## 
## REML criterion at convergence: 168.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.75223 -0.41953 -0.02529  0.44228  2.43853 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Bee.ID   (Intercept) 1.0672   1.0331  
##  Residual             0.3623   0.6019  
## Number of obs: 60, groups:  Bee.ID, 30
## 
## Fixed effects:
##                                   Estimate Std. Error t value
## (Intercept)                        2.34652    1.52140   1.542
## stroke.amplitude_scaled            1.12885    0.73555   1.535
## wbf.aud._scaled                    0.82929    0.21120   3.927
## percLoad                           5.34378    1.00541   5.315
## Treatment.order2                  -0.02245    0.18159  -0.124
## size_pca1_scaled                   2.79423    0.26035  10.732
## stroke.amplitude_scaled:percLoad  -1.08484    0.52736  -2.057
## Treatment.order2:size_pca1_scaled -0.44257    0.16981  -2.606
## 
## Correlation of Fixed Effects:
##             (Intr) strk._ wbf.._ percLd Trtm.2 sz_p1_ st._:L
## strk.mpltd_ -0.217                                          
## wbf.d._scld  0.431 -0.517                                   
## percLoad    -0.986  0.297 -0.485                            
## Trtmnt.rdr2  0.231 -0.186  0.500 -0.296                     
## sz_pc1_scld -0.027 -0.233  0.468 -0.002  0.218              
## strk.mpl_:L  0.514 -0.941  0.573 -0.587  0.251  0.182       
## Trtmn.2:_1_  0.281 -0.222  0.038 -0.291  0.016 -0.339  0.287
# write output
summary(mm11)$coefficients  
##                                      Estimate Std. Error    t value
## (Intercept)                        2.34652139  1.5213993  1.5423442
## stroke.amplitude_scaled            1.12884653  0.7355504  1.5346963
## wbf.aud._scaled                    0.82929053  0.2111955  3.9266482
## percLoad                           5.34378104  1.0054092  5.3150309
## Treatment.order2                  -0.02245326  0.1815899 -0.1236482
## size_pca1_scaled                   2.79422560  0.2603545 10.7323894
## stroke.amplitude_scaled:percLoad  -1.08484103  0.5273643 -2.0570998
## Treatment.order2:size_pca1_scaled -0.44257485  0.1698069 -2.6063418
write.csv( summary(mm11)$coefficients, file = "resp_oth_Coefs_percLoad.csv" )

# diagnostics
# qq plot
qqnorm(resid(mm11), main = "")
qqline(resid(mm11)) # good, but two outliers

# residual plot
plot(fitted(mm11), resid(mm11), xlab = "fitted", ylab = "residuals")
abline(0,0)

# QQPlot for group-level effects
qqnorm(ranef(mm11)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(mm11)$Bee.ID[[1]]) # looks good

This model is fundamentally different from the previous one. It’s saying, when we account for differences in wingbeat freq, stroke amplitude, and bee size, we sill see an effect of percent load on metabolic rate.

OR, we could say that holding all other variables constant, an increase in wingbeat frequency is associated with an increase in metabolic rate.

REFREF: Analyze for load (or percent load)

Analyze total mass is associated with metabolic rate ? (REFREF)